clear; close all; addpath('../utilities');


dc      = xlsread('../data/consumption.xlsx','H15:H91'); % 1947-2006 data
start   = 19; % 1947-48 growth rate
data    = [mean(dc(start:end)),std(dc(start:end)),skewness(dc(start:end)),kurtosis(dc(start:end)),corr(dc(start:end-1),dc(start+1:end))]';

p       = 0.878776;

mu1     = 0.003060;
mu2     = 0.000453;

sig1    = 0.001500;
sig2    = 0.002768;
sig3    = 0.007034;

rng('default')

DD           = 44; % mean disaster duration in month
i            = 1:100000;
duration_fun = @(p)sum((1-p)*p.^(i-1).*i);
p3          = fzero(@(p)duration_fun(p)-DD,0.99);

ED          = 0.834185141; % mean drop in consumption (simple growth rate)
i           = 1:100000;
mean_fun    = @(m3)sum((1-p3)*p3.^(i-1).*exp(i*m3+i*sig3^2/2));
mu3         = fzero(@(m)mean_fun(m)-ED,-0.0038);

pd          = fzero(@(x)solve_P(x,x,p,p,p3)-(44/12)/(2019-1929),1e-3);

mu_C        = [mu1;mu2;mu3];
sig_C       = [sig1;sig2;sig3];

P = [ p*(1-pd)      (1-p)*(1-pd)  pd;
      (1-p)*(1-pd)  p*(1-pd)      pd;
      (1-p3)/2      (1-p3)/2      p3];
  
stat        = P^1e6;
stat        = stat(1,:)';

P_nodis = [p,1-p;1-p,p];

Tsim       = (2006-1946)*12;
Psim       = 10000;
gA         = NaN(Tsim/12-1,Psim);

for i=1:Psim
    s0         = randi(2);
    [~,states] = sim_markov_chain_fast(P_nodis,Tsim,1,s0,1:2); % impose no disaster
    epsC       = randn(Tsim,1);
    gC         = mu_C(states) + sig_C(states).*epsC;
    C          = exp(cumsum(gC));
    C          = reshape(C,[12,Tsim/12]);
    Ca         = sum(C)';
    gA(:,i)    = log(Ca(2:end)./Ca(1:end-1));
end
gA = gA(:);

vol_gA = std(gA);
ac1_gA = (mean(gA(1:end-1,:).*gA(2:end,:)) - mean(gA(1:end-1,:)).*mean(gA(2:end,:)))./vol_gA.^2;

model = [mean(gA) vol_gA skewness(gA) kurtosis(gA) ac1_gA]';

table = [mu_C'*12; sig_C'*sqrt(12); [p, p3, pd*12]];
col_labels = {'$\mu_C \times 12$','$\sigma_C \times \sqrt{12}$','$P$'};
row_labels = {'1','2','3'};
Opts.Style = 'jf'; Opts.digit='r'; format = '%.4f';
disp(Mat2TexC(table,format,col_labels,row_labels,[],Opts,'Panel A: Parameters'))

table = [DD; (1-ED)*100; stat(3)*100; sig_C(3)*sqrt(12)*100];
col_labels = {'Average duration (months)','Average cumulative consumption drop','Unconditional probability','Conditional volatility'};
row_labels = {};
Opts.Style = 'jf'; Opts.digit='r'; format = '%.2f';
disp(Mat2TexC(table,format,col_labels,row_labels,[],Opts,'Panel B: Crises Moments'))

table = [data model];
col_labels = {'Mean','Std. Dev.','Skewness','Kurtosis','Auto-Corr.'};
row_labels = {'Data','Model'};
Opts.Style = 'jf'; Opts.digit='r'; format = '%.4f';
disp(Mat2TexC(table,format,col_labels,row_labels,[],Opts,'Panel C: Post-War Moments'))

